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INTRODUCTION 

The  I  low  around  bends  of  various  cross  section  and  through  other  curved  flow  pa. 
sages  is  of  considerable  importance  in  internal  flow  applications.  A  very  common 
example  is  that,  of  flow  through  an  elbow  or  other  section  of  curved  pipe  connecting  tw  ■ 
sections  of  straight  pipe.  Other  examples  include  that  of  flow  in  curved  ducts  of  re< 
(angular  or  other  cross  section  and  flow  in  more  general  smoothly-curved  passages. 
Although  such  geometries  are  unavoidable  in  many  instances  to  achieve  compact  and/or 
lightweight  designs,  they  are  sources  of  complex  secondary  flows,  losses,  and  variable 
heal  transfer  rates  which  may  affect  overall  performance  or  present  other  design  con¬ 
straints.  In  other  instances,  the  presence  of  geometries  which  produce  high  total, 
pressure  losses  or  heat  transfer  may  be  quite  intentional,  but  accompanied  by  a  need 
lor  deLailcd  understanding  and/or  predictive  techniques. 

The  underlying  physicaL  mechanisms  present  in  flows  of  this  type  are  clearly 
elucidated  by  secondary  flow  theory  (reviewed  for  example  by  Hawthorne  [1,2]  and  Hor1o<  l 
&  Lakshin  1  narayana  [3]).  In  its  simplest  form,  the  secondary  flow  is  generated  by  tinn¬ 
ing  a  primary  flow  in  which  viscous  or  other  forces  upstream  have  produced  a  non-zero 
velocity  gradient  normal  to  the  plane  of  curvature.  Fluid  with  above  (/below)  average 

momentum  migrates  to  the  outside  (/inside)  of  the  bend  as  a  result  of  the  radial  pres¬ 

sure  gradients  produced  by  turning  the  flow.  This  phenomenon  is  quantified  by  secondary 
Mow  theory  as  the  generation  of  streamwise  vorticity  from  transverse  vorticity  which 
lias  been  produced  upstream.  The  secondary  velocities  are  usually  quite  large  (a  sub- 

si.  ml ia I  fraction  of  the  primary  velocity)  and  are  influenced  by  several  factors  such 

as  I  low  deflection,  strength  of  inlet  vorticity,  and  thickness  of  shear  regions. 

Although  the  secondary  flow  is  generated  by  an  inviscid  mechanism,  its  strength 
and  subsequent  development  are  influenced  in  varying  degrees  by  viscous  effecls. 
Furthermore,  since  strong  deflections  may  occur  over  a  short  distance,  such  flows  arc 
usually  of  a  transition  type  and  never  become  fully  developed  or  assume  any  convenient 
similarity  form.  By  analogy  with  external  flows,  the  flow  often  behaves  an  an  inviscid 
flow  in  a  central  core  region,  with  viscous  effects  limited  to  regions  near  solid 
boundaries.  Unlike  most  external  flow  problems,  however,  the  inviscid  core  region  is 
often  not  an  irrotational  potential  flow  but  is  rotational  with  a  complex  interaction 
occurring  between  the  two  flow  regions.  Furthermore,  as  the  flow  passes  through  suc¬ 
cessive  flow  passages,  new  viscous  and  thermal  boundary  layers  develop  beneath  pre¬ 
vious  boundary  layers,  and  the  distinction  between  rotational  inviscid  and  viscous 
boundary  layer  regions  becomes  tenuous. 
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Previous  Work 


A  number  of  previous  investigations  have  computed  flow  through  strongly  curved 
ducts  using  approximate  forms  of  the  governing  equations. 

Pratap  and  Spalding  [A]  have  considered  turbulent  flow  in  a  strongly  cut  veil  duct 
using  their  "partially  parabolic"  calculation  procedure  and  a  two-equafcion/wai l -function 
turbulence  model.  Ghia  and  Sokhey  [5]  have  computed  laminar  flow  in  ducts  of  strong 
curvature  using  a  parabolized  form  of  the  Navier-Stokes  equations.  Kreskovsky,  Briley 
and  McDonald  [6]  have  recently  employed  an  approximate  initial-value  analysis  for 
viscous  primary  and  secondary  flows  to  compute  laminar  and  turbulent  flow  in  strongly 
curved  ducts,  using  a  one-equation  turbulence  model  with  viscous  sublayer  resolution. 
Humphrey,  Taylor  and  Whitelaw  [7,8]  have  obtained  laser-Doppler  anemometry  measurements 
of  laminar  and  turbulent  flow  in  a  square  duct  of  strong  curvature,  and  also  performed 
numerical  calculations  using  a  version  of  the  ful ly-elliptic  calculation  procedure 
developed  at  Imperial  College  by  Gosman,  Pun,  Patankar  and  Spalding.  Further  extensive 
calculations  including  heat  transfer  effects  have  been  made  recently  by  Yee  and 
Humphrey  [9]. 

In  view  of  the  general  complexity  of  flow  in  curved  ducts  and  pipes,  including 
the  likely  occurrence  of  flow  separation,  analysis  by  numerical  solution  of  the 
Navier-Stokes  equations  is  both  attractive  and  viable  provided  efficient  numerical 
methods  are  used.  This  approach  was  recently  employed  by  the  authors  [10,11]  co  study 
laminar  and  turbulent  flow  in  curved  ducts,  pipes,  and  two-dimensional  channels.  This 
study  concentrated  on  the  square-duct  geometry  and  flow  conditions  for  which  detailed 
measurements  have  been  obtained  by  Taylor,  Whitelaw  and  Yianneskis  [12]  .  In  this  study, 
mesh  resolution  tests  and  validation  of  inf low/outf low  boundary  conditions  were  per- 
l'o nned  on  the  two-dimensional  channel  problem.  Solutions  for  laminar  and  turbulent 
flow  on  corresponding  square  ducts  having  the  same  curvature  as  the  channels  were 
computed  and  compared  with  measurements  from  [12].  Good  qualitative  and  reasonable 
quantitative  agreement  between  solution  and  experiment  was  obtained. 

In  the  present  investigation,  further  extensive  calculations  were  performed 
using  the  same  computational  method  as  in  [10,11].  Six  flow  cases  including  both 
laminar  and  turbulent  flow  in  curved  pipes  and  square  ducts  are  evaluated  and  compared 
with  available  experimental  measurements.  The  laminar  square  duct  of  curvature  ratio 
2.3  (radius  of  curvature  R  to  hydraulic  diameter  d)  considered  in  [10,11]  was  recomputed 
to  test  its  sensitivity  to  the  grid  distribution  and  other  solution  parameters.  Two 
turbulent  calculations  for  this  s.tir._  geometry,  but  having  different  inlet  boundary  layer 
thicknesses  are  evaluated.  One  of  these  cases  has  a  fu, ly-developed  inflow  and  was 
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chosen  because  it  is  one  of  the  evaluation  cases  for  the  1980-81  Stanford  Conference 
on  Complex  Turbulent  Flows  for  comparison  of  computation  and  experiment.  Two  pipe 
bends  are  considered,  one  laminar  and  one  turbulent,  each  having  a  curvature  ratio  of 
2.8.  These  are  compared  with  the  recent  measurements  of  Taylor,  Whitelaw  and  Yianneskis 
[13].  Finally,  turbulent  flow  in  a  oharp  elbow  of  curvature  ratio  1.0  is  considered. 
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THE  PRESENT  APPROACH 


The  present  approach  is  the  numerical  solution  of  the  compressible  Reynolds- 
averaged  Navier-Stokes  equations  in  the  low  Mach  number  regime  (M  =  0.05)  for  which 
they  approximate  the  flow  of  a  liquid.  The  governing  equations  are  solved  using  an 
efficient  and  noniterative  time-dependent  linearized  block  implicit  (LBI)  scheme.  With 
proper  treatment  of  boundary  conditions,  this  algorithm  provides  rapid  convergence  which 
is  not  significantly  degraded  by  the  extreme  local  mesh  resolution  which  is  believed  to 
be  necessary  for  the  near-wall  sublayer  region  in  turbulent  flows.  The  computational 
method  has  been  described  by  the  authors  in  a  previous  study  of  flow  in  duct  and  pipe 
bends  [10,11],  For  completeness,  an  updated  account  of  the  method  is  included  here. 

Governing  Equations  and  Coordinate  System 

The  compressible  Navier-Stokes  equations  in  general  orthogonal  coordinates  are 
solved  using  analytical  coordinate  data  for  a  system  of  coordinates  aligned  with  the 
duct  geometry.  The  compressible  time-dependent  Navier-Stokes  equations  are  written  in 
orthogonal  coordinates  in  the  form  given  by  Hughes  and  Gaylord  [14],  Terms  which  are 
identically  zero  in  the  coordinate  system  used  are  omitted.  The  first-derivative  flux 
terms  are  written  in  conservation  form,  and  for  economy  the  stagnation  enthalpy  is 
assumed  constant.  The  definition  of  stagnation  enthalpy  and  the  equation  of  state 
for  a  perfect  gas  can  then  be  used  to  eliminate  pressure  and  temperature  as  dependent 
variables,  and  solution  of  the  energy  equation  is  unnecessary.  The  continuity  and 
three  momentum  equations  are  solved  with  density  and  the  three  velocity  components 
aligned  with  the  coordinates  as  dependent  variables. 

The  coordinate  system  is  shown  in  Fig.  1  and  consists  of  an  axial  coordinate 
parallel  to  a  curved  duct  centerline  (which  lies  in  the  Cartesian  x-y  plane),  and 
general  orthogonal  coordinates  x2,  x3  in  transverse  planes  normal  to  the  centerline. 
Straight  extensions  are  included  upstream  and  downstream  of  the  bend.  If  the  axial 
coordinate  x^  denotes  distance  along  the  centerline  and  if  K(x^)  =  1/R(x^)  is  the 
centerline  curvature,  then  the  metric  scale  factor  h^  for  the  axial  coordinate  direction 
is  given  by  h^  ■  1  +  K(x^)  AR(x2>  x^) ,  where  AR  £  r-R  is  independent  of  x^.  The  trans¬ 
verse  metric  factors  are  given  by  h2  ■  h3  ■  1  for  rectangular  (Cartesian)  cross  sections 
and  by  h2  ■  1,  hj  *  x2  for  circular  (polar)  cross  sections.  The  quantity  AR  is  given 
by  AR  »  x3  and  AR  -  x2  cos  x3  for  Cartesian  and  polar  cross  sections,  respectively. 

The  centerline  curvature  K  is  discontinuous  when  a  straight  segment  of  a  duct  joins  a 
constant  radius  segment.  To  remove  the  associated  coordinate  singularity,  the  flow 
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geometry  is  smoothed  over  an  axial  distance  of  o.  duct  diameter.  This  is  accomplished 
using  a  fifth-order  polynomial  variation  in  K  which  matches  function  value,  first  and 
second  derivatives  of  K  at  the  end  points  of  the  smoothing  region.  Analytical  coordinate 
transformations  due  to  Roberts  [15]  were  used  to  redistribute  grid  points  and  thus 
improve  resolution  in  shear  layers.  Derivatives  of  geometric  data  were  determined 
analytically  for  use  in  the  difference  equations.  A  nonunifora  mesh  distribution  was 
employed  for  the  axial  coordinate  to  concentrate  grid  points  inside  the  curved 
section  of  the  duct  or  pipe.  Representative  grids  for  the  cross  sectional  planes  are 
shown  in  Fig.  2. 


Turbulence  Model 


The  turbulence  model  used  falls  into  the  category  of  one-equation  turbulence 
models  discussed  by  Launder  and  Spalding  [16],  and  parallels  the  method  given  by 
Shararoth  and  Gibeling  [17].  This  model  requires  solution  of  a  single  partial  dif¬ 
ferential  equation  governing  turbulence  kinetic  energy  k,  in  conjunction  with  a 

specified  length  scale  A.  A  turbulent  viscosity  p  is  obtained  from  the  Prandtl- 

1/2  C 

Kolmogorov  constitutive  relationship  pt  ot  £  k  .  The  turbulent  viscosity  pfc  is 
assumed  to  be  isotropic,  and  the  stress  tensor  in  the  ensemble-averaged  equation  is 
determined  by  adding  the  turbulent  viscosity  to  the  molecular  viscosity  p  to  obtain 
the  total  effective  viscosity  pg  ■  p  +  y^.  Given  an  estimate  of  the  length  scale  at 
the  edge  of  the  viscous  layer  and  its  streamwise  growth  rate,  a  distribution  of  mixing 
length  within  the  viscous  layer  can  be  obtained  from  semi-empirical  relationships 
widely  used  in  two-  and  three-dimensional  boundary  layer  calculations.  To  estimate 
the  axial  growth  rate  of  the  shear  layer,  a  boundary  layer  momentum  integral  procedure 
which  neglects  axial  curvature  but  includes  blockage  effects  is  used.  The  final 
turbulence  model  provides  for  resolution  of  the  viscous  sublayer  region  near  walls. 

The  equation  governing  turbulence  kinetic  energy  is  given  by  Launder  and 
Spalding  [18]  for  Cartesian  coordinates  and  is  rewritten  in  the  present  orthogonal 
coordinate  system.  The  turbulent  viscosity  is  obtained  from  the  Prandtl-Kolmogorov 
relation 


pk  /e 


„  1/A  1/2. 

pk  £ 


(1) 


where  the  dissipation  rate  e  is  given  by 


e 


(2) 
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For  high  Reynolds  number  flow,  is  approximately  0.09.  For  low  Reynolds  number  and 
in  the  viscous  sublayer,  is  given  a  prescribed  dependence  on  turbulence  Reynolds 
number  R  in  the  manner  suggested  by  McDonald  and  Fish  [19]  for  transitional  boundary 
layer  flows.  In  [19],  the  turbulence  Reynolds  number  was  appropriately  defined  as 
an  integral  average  across  the  boundary  layer.  Here  and  following  [17],  it  If  convenient 
to  define  Rt  as  the  local  ratio  of  turbulent  to  laminar  viscosity  R^  =  ut/y.  The 
structural  coefficient  c^  is  given  in  [19]  as  c^  ■  4(a^)2,  where 

a!  =  a0  f(V  /  (1  +  6,66  a0(f(RT)_1)1  (3> 


where  a  *  0.0115  and 
o 

f(RT>  - 


0.22 


6.81  Rt 


+  6.143 


R  *  1 

T 

R  -  40 


(4) 


with  a  cubic  polynomial  curve  fit  for  values  of  R^  between  1  and  40. 

It  remains  to  specify  a  length  scale  distribution  appropriate  for  the  problem 
under  consideration  The  length  scale  distribution  is  adapted  from  previous  turbulence 
models  for  turbulent  boundary  layers  and  taken  to  be  the  conventionally  defined  mixing 
length  of  Prandtl.  The  distribution  of  mixing  length  given  by  McDonald  and  Fish  [19] 
has  proven  effective  for  a  wide  range  of  two-dimensional  turbulent  boundary  layers  and 
is  easily  adapted  for  present  use.  This  distribution  is  given  by 


£  =  ^  £  tanh  (xd/£  )  (5) 

OO  CJO 

where  £  is  mixing  length,  im  is  an  outer-region  value  of  mixing  length,  d  is  distance 
from  the  wall,  k  is  the  von  Karman  constant  (taken  here  as  0.41),  and  is  a 
sublayer  damping  function  given  by 

2>=  P1/2  [(d+-23)/8J  (6) 

Here,  P  is  the  normal  probability  function  and  d+  is  defined  by  d+  =  d(r  /^)1/2/v, 
where  t  is  shear  stress  and  v  is  kinematic  viscosity.  For  equilibrium  boundary  layers, 
is  observed  to  have  a  constant  value  of  about  0.09  6,  where  6  is  the  local  boundary 
layer  thickness. 

The  length  scale  distribution  of  Eq.  (5)  is  adapted  for  present  use  by  taking  d  as 
distance  to  the  nearest  wall  and  by  assigning  £ a  distribution  based  on  two-dimensional 
momentum  integral  estimates  of  the  boundary  layer  growth  rates.  The  computed  estimates 
of  boundary  layer  growth  were  obtained  from  a  simple  integral  prediction  scheme  which 
neglects  axial  curvature,  rather  than  an  attempt  to  scan  the  intermediate  transient 
solutions  of  the  Navier-Stokes  equations  to  determine  some  ill-defined  boundary  layer 
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thickness.  Assuming  a  1/7  power  law  velocity  profile  and  a  skin  friction  law 

1/4 

c^/2  ■  0.0225  (v/u  5)  ,  the  momentum  integral  equation  for  a  two-dimensional 

r  ^ 

straight  channel  can  be  written  as 


.  dfi  .  _6_  du« 

hl  dx  +  u  dx 
e 


Of/2 


(7) 


where  for  the  1/7  power  profile:  h.  =  7/72  and  h2  =  22/72.  An  analogous  equation  is 

easily  derived  for  flow  in  a  straight  pipe.  The  free  stream  velocity  is  then  related 

to  the  be  'r  1ary  layer  thickness  by  an  ancillary  formula  which  approximates  the  blockage 

effects  a  .  .jera  ed  with  internal  flow.  This  relationship  between  6  and  u  is  inserted 

£ 

intf  .  x/)  prior  to  solution,  ilq.  (7)  is  solved  using  a  second-order  linearized 
Jj  z  difference  scheme  described  in  [23].  The  outer  or  maximum  mixing  length  scale 
a-  'ci  ted  with  c.  ~h  axial  location  is  determined  from  the  formula  Zm  =  0.09  <$,  where 
<i  is  '  »  c  imputed  value  from  Eq.  (7). 


Physical  Boundary  and  Initial  Conditions 

Toe  computational  domain  includes  the  curved  portion  of  the  duct  or  pipe  compris¬ 
ing  the  90-degree  bend  and  short  straight  extensions  upstream  and  downstream  of  the 
bend.  The  computational  domain  is  thus  embedded  within  a  larger  overall  flow  system 
and  requires  inflow  and  outflow  boundary  conditions  which  adequately  model  the  inter¬ 
face  between  the  computed  flow  and  the  remainder  of  the  flow  system.  The  inflow/ 
outflow  conditions  used  are  derived  from  an  assumed  flow  structure  and  are  chosen  to 
provide  inflow  with  prescribed  stagnation  pressure  (and  stagnation  enthalpy)  in  an 
inviscld  core  region  and  with  a  given  axial  velcoity  profile  shape  in  the  inflow 
boundary  layer,  and  to  provide  outflow  with  a  prescribed  distribution  of  static  pres¬ 
sure  in  the  cross  section.  These  boundary  conditions  are  compatible  both  with  an 
inviscid  characteristics  analysis  and  with  the  physical  process  by  which  most  such 
flows  are  established.  Physically,  a  duct  flow  is  often  established  by  supplying 
air  of  a  given  stagnation  pressure  and  temperature  and  exhausting  the  duct  at  a  given 
static  pressure.  The  mass  flux  through  the  duct  may  then  vary  with  time  until  a  steady 
state  is  achieved,  at  which  the  mass  flux  is  determined  as  a  balance  between  these 
inflow/outflow  quantities  and  viscous  and  thermal  effects  within  the  duct.  By  choosing 
stagnation  pressure  and  temperature  at  inflow  and  static  pressure  at  outflow  as  the 
dominant  boundary  conditions,  the  present  solution  procedure  allows  both  velocity  and 
static  pressure  to  vary  with  time  at  the  inflow  boundary,  C3  is  consistent  with  this 


physical  process.  As  a  consequence,  pressure  waves  are  transmitted  upstream  through 
the  inflow  boundary  during  the  transient  flow  process  and  are  not  reflected  back  into 
the  computational  domain.  The  reflection  of  pressi.re  waves  at  an  inflow  boundary 
when  velocity  and  pressure  are  fixed  in  time  has  often  been  cited  as  a  cause  of 
either  instability  or  slow  convergence  in  other  investigations.  These  boundary 
conditions  are  discussed  in  more  detail  in  [20].  The  specific  treatment  of  initial 
and  boundary  conditions  used  here  is  outlined  below. 

The  initial  and  boundary  conditions  are  devised  from  estimates  of  the  potential 
flow  velocity  U^(x^,  x^,  x^)  for  the  duct,  a  mean  boundary  layer  thickness  <$(x^)  for 
shear  layers  on  transverse  duct  or  pipe  walls,  and  finally  from  an  estimate  of  the 
blockage  correction  factor  B(x^)  for  the  core  flow  velocity  due  to  the  boundary  layer 
growth.  The  potential  flow  velocity  is  approximated  as  uniform  flow  in  straight 

segments  of  the  duct  and  as  Cr  ^  in  curved  segments.  Taking  C  as  (R  -  R.  )/£n(R  /R. ) 

2  221/2  olol. 

and  d  / 8 [R—  (R  -d  /4)  ]  leads  to  a  unit  mass  flux  for  rectangular  and  circular  cross 

sections,  respectively.  Here,  R^  and  Rq  are  the  radii  to  the  inner  and  outer  walls 

of  the  duct,  R  =  (R^  +  RQ)/”«  and  d  is  the  pipe  diameter  (cf.  Fig.  1).  Distributions 

of  6(x^)  and  Bfx^)  are  determined  by  recourse  to  the  simple  one-dimensional  momentum 

integral  analysis  mentioned  in  describing  the  turbulence  model.  Finally,  a  shear 

layer  velocity  profile  shap«  f(y/6),  o-f-1  is  chosen  for  each  problem,  where  y  is  a 

parameter  indicative  of  distance  from  a  wall.  For  laminar  flow  cases,  von  Karman 

polynomial  velocity  profiles  were  used.  For  turbulent  flow  cases,  these  velocity 

profile  shapes  were  taken  from  the  analytical  fit  of  Musker  [21]  to  the  Coles  type 

of  profile  which  matches  <5  and  c^  from  the  momentum- integral  calculations.  The 

initial  values  of  velocity  components  u^,  u2,  are  given  by 

■  Uj  B(x1)  f  [y/6(x1>] 

u2  =  u3  =  0 

In  the  pipe  flow  calculations,  estimates  of  the  inflow  radial  velocity  u2  were 
available  and  were  used  as  initial  conditions.  The  details  of  this  procedure  are 
not  critical  and  are  omitted,  since  except  for  tho  shape  and  thickness  of  the  inlet 
boundary  layer  profile,  these  results  serve  only  as  a  convenient  method  for  selecting 
initial  conditions.  A  reasonably  accurate  estimate  for  the  pressure  drop  which  will 
produce  the  desired  flow  rate  must  be  made  using  any  convenient  source,  such  as  a 
Moody  diagram,  data  correlations,  momentum  integral  analysis,  or  other  computed 
results.  A  smooth  axial  distribution  of  pressure  which  matches  this  pressure  drop 
is  then  assigned  aad  adjusted  to  approximate  local  curvature  of  the  flow  geometry. 
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This  completes  specification  of  the  initial  conditions.  It  is  noted  that  although 
these  initial  conditions  do  take  into  account  several  relevant  features  of  the  flow, 
the  important  effects  of  strong  secondary  flows  and  their  distortion  of  the  primary 
flow  are  completely  neglected.  The  resulting  initial  flow  is  thus  a  simple  but 
relatively  crude  approximation  to  the  final  flow  field. 

At  the  inflow  boundary,  a  "two-layer"  boundary  condition  is  devised  such  that 
stagnation  pressure  is  fixed  in  the  core  flow  region  and  an  axial  velocity  pro¬ 
file  shape  ■  f(y/6)  is  set  within  shear  layers.  Here,  ug  is  the  local  edge 

velocity  which  varies  with  time  and  is  adjusted  after  each  time  step  to  the  value 
consistent  with  pQ  and  the  local  edge  static  pressure  determined  as  part  of  the 
solution.  The  transverse  velocity  components  U2  and  u^  were  also  specified  at 
intlow.  For  the  pipe  flow  cases,  was  set  to  zero,  and  the  radial  velocity  Ug 
was  given  a  smooth  distribution  consistent  with  the  boundary  layer  thickness  and 
estimated  axial  pressure  drop.  For  the  duct  flow  cases,  and  u^  were  extrapolated 
until  the  initial  impulsive  transients  had  passed  and  were  then  held  fixed.  The 
remaining  inflow  condition  is  3  c^/3n  =  g  (x^,  x^) ,  where  n  denotes  the  normal 

coordinate  direction  and  c  is  pressure  coefficient  referred  to  reference  conditions. 

"  2 
The  quantity  g  is  computed  from  the  initial  conditions  with  c^  defined  as  l-(BUj)  , 

its  value  from  the  potential  flow  corrected  for  estimated  blockage.  For  outflow 

conditions,  the  static  pressure  is  impcjed,  and  second  normal  derivatives  cf  each 

velocity  component  are  set  to  zero.  At  no-slip  walls,  all  velocity  components  are 

set  to  zero,  and  the  remaining  condition  used  is  3p/3n  =  0,  where  p  is  pressure. 

The  condition  3p/3n  =  0  at  a  no-slip  surface  approximates  the  normal  momentum 

equation  to  order  Re  ^  foi  viscous  flow  at  high  Reynolds  number.  Finally,  the  three- 

dimensional  flow  cases  are  assumed  to  be  symmetric  about  the  plane  containing  the 

curved  duct  centerline,  and  symmetry  conditions  are  imposed  on  this  boundary. 


SOLUTION  PROCEDURE 
Background 


The  solution  procedure  employs  a  consistently-split  linearized  block  implicit 
(LBI)  algorithm  which  has  been  discussed  in  detail  by  the  authors  in  [22,23].  There 
are  two  important  elements  of  this  method: 

(1)  the  use  of  a  noniterative  formal  time  linearization  to  produce 
a  fully-coupled  linear  multidimensional  scheme  which  is  written 
in  "block  implicit"  form;  and 

(2)  solution  of  this  linearized  coupled  scheme  using  a  consistent 
"splitting"  (ADI  scheme)  patterned  after  the  Douglas-Gunn  [24] 

(1964)  treatment  of  scalar  ADI  schemes. 

The  method  is  thus  referred  to  as  a  split  linearized  block  implicit  (LBI)  shcerae. 
The  method  has  several  attributes: 

(1)  the  noniterative  linearization  is  efficient; 

(2)  the  fullv-coupled  linearized  algorithm  eliminates  instabilities 
and/or  extremely  slow  convergence  rates  often  attributed  to 
methods  which  employ  ad  hoc  decoupling  and  linearization 
assumptions  to  identify  nonlinear  coefficients  which  are  then 
treated  by  lag  and  update  techniques; 

(3)  the  splitting  or  ADI  technique  produces  an  efficient  algorithm 
which  is  stable  for  large  time  steps  and  also  provides  a  means 
for  convergence  acceleration  for  further  efficiency  in  computing 
steady  solutions; 

(4)  intermediate  steps  of  the  splitting  are  consistent  with  the 
governing  equations,  and  this  means  that  the  "physical"  boundary 
conditions  can  be  used  for  the  intermediate  solutions.  Other 
splittings  which  are  inconsistent  can  have  severe  difficulties  ii. 
satisfying  physical  boundary  conditions  [23]. 

(5)  the  convergence  rate  and  over?' '  efficiency  of  the  algorithm  are 
much  less  sensitive  to  mesh  refinement  and  redistribution  than 
algorithms  based  on  explicit  schemes  or  which  employ  ad  hoc 
decoupling  and  linearization  assumptions.  This  is  important  for 
accuracy  and  for  computing  turbulent  flows  with  viscous  sublayer 
resolution;  and 
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(6)  the  method  is  general  and  is  specifically  designed  for  the 
complex  systems  of  equations  which  govern  multiscale  viscous 
flow  in  complicated  geometries. 

This  same  algorithm  was  later  considered  by  Beam  and  Warming  [25],  but  the  ADI  split¬ 
ting  was  derived  by  approximate  factorization  instead  of  the  Douglas-Gunn  procedure. 
They  refer  to  the  algorithm  as  a  "delta  form"  approximate  factorization  scheme. 

This  scheme  replaced  an  earlier  non-delta  form  scheme  [26],  which  has  inconsistent 
intermediate  steps. 


Spatial  Differencing  and  Artificial  Dissipation 


The  spatial  differencing  procedures  used  are  a  straightforward  adaption  of 
those  used  in  [22]  and  elsewhere.  Three-point  central  difference  formulas  are  used 
for  spatial  derivatives,  including  the  first-derivative  convection  and  pressure 
gradient  terms.  This  has  an  advantage  over  one-sided  formulas  in  flow  calculations 
subject  to  "two-point"  boundary  conditions  (virtually  all  viscous  or  subsonic  flows), 
in  that  all  boundary  conditions  enter  the  algorithm  implictly.  In  practical  flow 
calculations,  artificial  dissipation  is  usually  needed  and  is  added  to  control  high- 
frequency  numerical  oscillations  which  otherwise  occur  with  the  central-difference 
fcrmula. 

In  the  present  investigation,  artificial  (anisotropic)  dissipation  terms  of 
the  form 


Z 

i 


(8) 


are  added  to  the  right-hand  side  of  each  (k-th)  component  of  the  momentum  equation, 
where  for  each  coordinate  direction  x j ,  the  artificial  diffusivity  d^  is  positive 
and  is  chosen  as  the  larger  of  zero  and  the  local  quantity  y  (a  Re^x  -1)/Re.  Here, 
the  local  cell  Reynolds  number  Re.  for  the  j-th  direction  is  defined  by 

j 


Re 


Ax. 


Re  | pUj |  AXj/ye 


(9) 


This  treatment  lowers  the  formal  accuracy  to  0  (Ax) ,  but  the  functional  form  is  such 
that  accuracy  in  representing  physical  shear  stresses  in  thin  shear  layers  with  small 
normal  velocity  is  not  seriously  degraded.  This  latter  property  follows  from  the 
anisotropic  form  of  the  dissipation  and  the  combination  of  both  small  normal  velocity 
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and  small  grid  spacing  in  thin  shear  layers.  A  value  of  0.5  was  used  for  a  in  the 
present  calculations.  Values  lower  than  0.5  have  been  used  to  good  effect  in  two 
space  dimensions  [27,28],  but  it  has  not  yet  been  possible  to  investigate  the  role 
of  smaller  a  values  in  three  space  dimensions,  although  it  is  believed  that  lower 
values  of  o  would  also  be  beneficial  here. 

Split  LBI  Algorithm 
Linearization  and  Time  Differencing 

The  system  of  governing  equations  to  be  solved  consists  of  five  equations: 
continuity,  three  components  of  momentum  and  the  turbulence  kinetic  energy  equation 
in  five  dependent  variables:  p,  u^,  Uj,  u^,  and  k.  Using  notation  similar  to  that 
in  [22],  at  a  single  grid  point  this  system  of  equations  can  be  written  in  the 
following  form: 

3  H($)/3t  -  D<$)  +  S($)  (1Q) 

where  <t>  is  the  column-vector  of  dependent  variables,  H  and  S  are  column-vector 
algebraic  functions  of  $,  and  D  is  a  column  vector  whose  elements  are  the  snatial 
differential  operators  which  generate  all  spatial  derivatives  appearing  in  the 
governing  equation  associated  with  that  element. 

The  solution  procedure  is  based  on  the  following  two-level  implicit  time- 
difference  approximations  of  (10) : 

(Hn+1  -  Hn)/At  -  B(Dn+1  +  Sn+1)  +  (1-B)  (Dn  +  Sn)  (11) 

where,  for  example,  Hn+^  denotes  H ( 4>n~*"^ )  and  At  ■  tn+^  -Tn.  The  parameter  8 
(0.5  -  0  -  1)  permits  a  variable  time- centering  of  the  scheme,  with  a  truncation 
error  of  order  [At2,  (0  -  1/2)  At]. 

A  local  time  linearization  (Taylor  expansion  about  $n)  of  requisite  formal 
accuracy  is  introduced,  and  this  serves  to  define  a  linear  differential  operator  L 


(Cf.  [22]) 

such  that 

Similarly, 

Dn+1  =  Dn  +  Ln  (<frn+1  -  <f>n)  +  0  (At2) 

(12) 

Hn+1  -  Hn  +  (3H/3<J>)n  (<frnrl  -  $")  +  0  (At2) 

(13) 

sn+l  e  sn  +  os/g^)0  ($n+1  _  vn)  +  0  (At2) 

(1A) 
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Eqs.  (12-14)  are  insetted  into  Eq.  (11)  to  obtain  the  following  system  which  is 
linear  in  4>n+'*' 

(A  -  BAt  ,n)  ($n+A  -  $n)  -  At  (Dn  +  Sn)  (15) 

and  which  is  termed  a  linearized  block  implicit  (LBI)  scheme.  Here,  A  denotes  a 
square  matrix  defined  by 

A  =  (3H/3$)n  -  BAt  (3S/3$)n  (16) 

Eq.  (15)  has  0  (At)  accuracy  unless  H  5  <£,  in  which  case  the  accuracy  is  the  same  as 
Eq.  (11). 

Special  Treatment  of  Diffusive  Terms 

The  time  differencing  of  diffusive  terms  is  modified  to  accomodate  cross-derivative 
terms  and  also  turbulent  viscosity  and  artificial  dissipation  coefficients  which  depend 
on  the  solution  variables.  Although  formal  linearization  of  the  convection  and  pres¬ 
sure  gradient  terms  and  the  resulting  implicit  coupling  of  variables  is  critical  to 
the  stability  and  rapid  convergence  of  the  algorithm,  this  does  not  appear  to  be 
important  for  the  turbulent  viscosity  and  artificial  dissipation  coefficients.  Since 

the  relationship  between  y  .ind  d,  and  the  mean  flow  variables  is  not  conveniently 

®  2  n 
linearized,  these  diffusive  coefficients  are  evaluated  explicitly  at  t  during  each 

time  step.  Notationally ,  this  is  equivalent  to  neglecting  terms  proportional  to 

3y  /3*  or  3d . / 3 <p  in  Ln,  which  are  formally  present  in  the  Taylor  expansion  (12),  but 

®  J  n  n 

retaining  all  terms  proportional  to  or  dj  in  both  L  and  D  . 

It  has  been  found  through  extensive  experience  that  this  has  little  if  any  effect 
on  the  performance  of  the  algorithm.  This  treatment  also  has  the  added  benefit  that 
the  turbulence  model  equatiors  (in  this  instance  the  turbulence  kinetic  energy  equation) 
can  be  decoupled  from  the  system  of  mean  flow  equations  by  an  appropriate  matrix  par¬ 
titioning  (cf.  [23])  and  solved  separately  in  each  step  of  the  ADI  solution  procedure. 
This  reduces  the  block  size  of  the  block  tridiagonal  systems  which  must  be  solved  in 
each  step  and  thus  reduces  the  computational  labor. 

In  addition,  the  viscous  terms  in  the  present  formulation  include  a  number  of 
spatial  cross-derivative  terms.  Although  it  is  possible  to  treat  cross-derivative 
terms  implicitly  within  the  ADI  treatment  which  follows,  it  is  not  at  all  convenient 
to  do  so,  and  consequently,  all  cross-derivative  terms  are  evaluated  explicitly  at  tn. 
For  a  scalar  model  equation  representing  combined  convection  and  diffusion,  it  has 
been  shown  by  Bema  and  Warming  [29]  that  the  explicit  treatment  of  cross-derivative 
terms  does  not  degrade  the  unconditional  stability  of  the  present  algorithm.  To 
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preserve  notational  simplicity,  it  is  understood  that  all  cross-derivative  terms 
appearing  in  L°  are  neglected  but  are  retained  in  Dn.  It  is  important  to  note  that 
neglecting  terms  in  Ln  has  no  effect  on  steady  solutions  of  Eq.  (15),  since 
=  o  and  thus  Fq.  (15)  reduces  to  the  steady  form  of  the  equations: 

Dn  +  Sn  =  0.  Aside  from  stability  considerations,  the  only  effect  of  neglecting 
terms  in  Ln  is  to  introduce  an  0  (At)  truncation  error. 

Consistent  Splitting  of  the  LBI  Scheme 

To  obtain  an  efficient  algorithm,  the  linearized  system  (15)  is  split  using 
ADI  techniques.  To  obtain  the  split  scheme,  the  multidimensional  operator  L  is 
rewritten  as  the  sum  of  three  "one-dimensional"  sub-operators  (i  ■  1,  2,  3)  each 
of  which  contains  all  terms  having  derivatives  with  respect  to  the  i-th  spatial 
coordinate.  The  split  form  of  Eq.  (15)  can  be  derived  either  as  in  [22,  23]  by 
following  the  procedure  described  by  Douglas  and  Gunn  ] 24 ]  in  their  generalization 
and  unification  of  scalar  AnT  schemes,  or  using  approximate  factorization  as  in  [25]. 
For  the  present  system  of  equations,  the  split  algorithm  is  given  by 

(A  -  SAtLj)  ($*  -  $n)  -  At  (Dn  +  Sn)  (17a) 

(A  -  0AtL")  ($**-  *")  «  A  ($*  -  *n)  (17b) 

(A  -  0AtL^)  ($n+1-  $n)  -  A  ($**-  <frn)  (17c) 

where  <f>  and  <f>  ure  consistent  intermediate  solutions  [22,  ] 2 ] .  If  spatial  deriva¬ 
tives  appearing  in  and  D  are  replaced  by  three-point  difference  formulas,  as 

indicated  previously,  then  each  step  in  Eqs.  (17a-c)  can  be  solved  by  a  block- 
cridiagonal  elimination. 

Combining  Eqs.  (17a-c)  gives 

(A  -  BAtL?)  A_1  (A  -  BAtl^)  A_1  (A  -  BAtL^)  (4>n+1  -  4>n) 

1  2  J  (18) 

-  At  (Dn  +  Sn) 

2 

which  approximates  the  unsplit  scheme  (15)  to  0  (At  ).  Since  the  intermediate  steps 

are  also  consistent  approximations  for  Eq.  (15N ,  physical  boundary  conditions  can  be 

used  for  <|>  and  4>  [22,  23].  Finally,  since  the  L.  are  homogeneous  operators,  it 

1  n+1  * 

follows  from  Eqs.  (17a-c)  that  steady  solutions  have  the  property  that  $  “  <t>  m 

^  ^  n 

•-  9  and  satisfy 

Dn  +  Sn  -  0  (19) 

The  steady  solution  thus  depends  only  on  the  spatial  difference  approximation  used 
for  (19)  and  does  not  depend  on  the  solution  algorithm  itself. 


COMPUTED  RESULTS 


Solutions  for  three-dimensional  laminar  and  turbulent  flow  in  90-degree  bends 
having  strong  curvature  and  both  circular  and  square  cross  sections  are  presented 
here.  To  assess  the  degree  of  accuracy  obtainable  in  this  type  of  three-dimension. u 
flow  calculation,  it  is  obviously  valuable  to  have  experimental  measurements  for 
comparison.  Fortunately,  laser-Doppler  anemometry  measurements  have  recently  become 
available  [7,  8,  12,  13]  for  several  laminar  and  turbulent  square  duct  and  pipe  bends;. 
Computed  results  for  six  different  flow  cases  are  included  here,  and  the  relevant  flow 
parameters  describing  these  cases  are  given  in  Table  I.  Each  of  the  three  bend 
geometries  considered  is  shown  to  scale  in  Fig.  3.  The  cases  include  both  laminar 
and  turbulent  flow  in  a  pipe  bend  of  curvature  2.8  (cases  1  and  2)  and  in  a  square 
duct  of  curvature  2.3  (cases  4  and  5).  Measurements  are  given  by  Taylor,  White  law  .in' 
Yianneskis  both  for  the  pipe  bend  [13]  and  for  the  curved  duct  [12]. 
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The  other  cases  include  a  pipe  elbow  configuration  (Case  3),  which  has  a  curvature 
of  1.0,  but  is  otherwise  identical  to  Case  2,  and  a  curved  duct  (Case  6)  which  is 
identical  to  Case  5  except  that  it  has  a  nominally  fully-developed  inflow  velocity 
profile.  The  latter  flow  case  was  measured  by  Humphrey,  Taylor  and  Whitelaw  [8] 
and  was  Included  as  one  of  the  evaluation  cases  for  the  1980-81  Stanford  conference 
on  Complex  Turbulent  Flows. 

A  vast  amount  of  information  is  obtained  from  the  computation  of  six  three- 
dimensional  flow  cases,  and  only  selected,  representative  results  can  be  included 
here.  The  present  approach  to  this  problem  is  as  follows:  First,  the  computed 
axial  velocity  is  compared  with  experimental  measurements  in  each  of  the  five  cases 
for  which  data  is  available.  Next,  each  of  the  three  pipe  bend  cases  (1  -  3)  is 
presented  in  a  sequence  of  plots  designed  to  show  the  developing  flow  structure. 

The  structure  of  square-duct  flows  was  considered  in  more  detail  in  [10,  11]. 

Finally,  a  comparison  is  made  of  the  two  turbulent  curved  duct  flows  which  differ 
only  in  the  inflow  boundary  layer  thickness. 

Convergence  Rate 

Before  proceeding  to  the  discussion  of  computed  flow  behavior,  an  indication  is 
given  of  the  degree  and  rate  of  convergence  obtained  in  these  calculations.  Turbulent 
flows,  particularly  in  three  dimensions,  are  especially  difficult  to  compute  because 
of  the  very  high  degree  of  mesh  redistribution  required  to  resolve  viscous  sublayer 
regions.  Since  peak  secondary  flow  velocities  often  occur  within  or  near  the  viscous 
sublayer,  resolution  of  the  viscous  sublayer  is  believed  necessary  to  obtain  an 
adequate  representation  of  the  overall  flow,  but  poses  a  computational  problem  for 
which  rapid  convergence  is  difficult  to  achieve.  The  complicated  flow  structure  in 
curved  ducts  and  pipes,  which  includes  strong  secondary  flows  and  their  resulting 
distortion  of  the  primary  flow,  also  contributes  to  the  difficulty  of  these  flows. 

The  present  computational  method  employs  variable  tine  steps  to  accelerate 
convergence,  as  discussed  in  [30].  Several  procedures  for  time  step  selection  are 
presently  under  evaluation,  but  typically  small  time  steps  are  used  Initially  during 
the  elimination  of  impulsive  transients,  followed  by  larger  time  steps  and  the  cyclic 
use  of  a  sequence  of  time  steps.  In  addition,  a  smaller  time  step  was  used  near  the 
leading  edge  than  elsewhere  in  the  flow  field.  Small  time  steps  are  effective  in 
reducing  high  (spatial)  frequency  error  components,  while  large  time  steps  are 
effective  in  reducing  low  frequency  errors. 
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An  indication  of  both  the  degree  and  rate  of  convergence  obtained  for  the  three- 
pipe  fl^w  calculations  (which  are  representative  cases)  is  shown  in  Fig.  A.  For 
comparison,  a  curve  representative  of  the  convergence  behavior  typically  obtained  for 
two-dimensional  turbulent  flow  cases  is  also  shown.  Although  the  present  convergence 
rate  is'  noticeably  slower  than  that  observed  using  this  same  algorithm  in  other 
calculations,  it  is  adequate  for  present  purposes.  The  present  calculations  were 
terminated  as  shown  in  Fig.  A.  for  reasons  of  economy,  since  it  appeared  that  the 
observed  changes  in  the  solution  at  this  point  were  of  minor  significance  and  would 
not  alter  any  of  the  conclusions  drawn  from  these  results. 

Comparison  With  Experiment 

Computed  contours  of  axial  velocity  are  compared  with  experimental  measurements 
in  Figs.  5-9.  To  aid  in  the  comparison,  contour  values  shorn  for  the  computed 
solutions  are  the  same  as  has  been  shown  for  the  measuremtns.  The  computed  velocities 
were  normalized  by  the  mean  (bulk)  velocity  as  determined  by  a  second-order  numerical 
integration  of  the  computed  solution.  This  calculation  indicated  that  global  mass 
conservation  was  satisfied  by  all  computed  solutions  to  within  one  per  cent  at  each 
axial  location.  In  comparing  these  results  with  the  measurements,  the  normalized 
mass  flux  (as  evidenced  by  peak  and  centerline  velocity)  did  not  appear  to  match  in 
some  cases.  To  aid  in  the  comparison,  the  computed  contours  shown  in  Figs.  5-9 
were  renormalized  to  provide  a  nominal  matching  of  centerline  velocity  at  the  first 
measured  station.  Downstream  stations  were  adjusted  by  this  same  factor.  This 
renormalization  led  to  adjustment  of  the  computed  velocities  by  the  following  amounts: 
Case  2  -  3.6%,  Case  A  -  7.A%,  Case  6  -  A%.  No  adjustment  was  required  for  Case  1, 
and  the  computed  mass  flux  was  not  checked  for  Case  5. 

All  of  the  flow  cases  considered  here  develop  strong  secondary  flows  as  a  result 
of  turning  and  these  secondary  flows  result  in  considerable  distortion  of  the  primary 
flow,  particularly  near  the  inside  of  the  bend.  This  distortion  of  the  primary  flow 
is  the  principal  method  of  judging  the  agreement  between  measured  and  computed  results 
in  Figs.  5-9.  Both  a  comparison  of  the  duct  flows  with  two-dimensional  channel  flow 
having  the  same  curvature  and  also  limited  mesh  refinement  studies  were  given  in  a 
previous  study  [10,  11 j . 

In  Fig.  5,  the  computed  results  for  laminar  pipe  flow  are  in  very  good  general 
agreement  with  the  measurements.  The  turbulent  pipe  flow  predictions  in  Fig.  6  also 
agree  reasonably  well  with  measurements,  although  the  predicted  flow  distortion  is 
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somewhat  less  localized  than  that  measured.  Since  the  laminar  calculations  employ 
no  physical  assumptions  beyond  those  of  the  Navier-Stokes  equations,  the  level  of 
agreement  is  indicative  of  the  numerical  truncation  error  and  experimental  error. 

The  difference  between  the  level  of  agreement  for  laminar  and  turbulent  predictions 
is  at  least  partly  an  indication  of  error  in  the  turbulence  model.  The  predictions 
for  the  three  duct  flow  cases  in  Figs.  7-9  are  in  good  qualitative  and  reasonable 
quantitative  agreement  with  the  measurements,  although  the  agreement  is  not  as  good 
as  was  obtained  in  the  pipe  flow  cases.  This  may  be  partly  due  to  the  increased 
difficulty  of  grid  resolution  (cf.  Fig.  2)  in  the  center  of  the  duct,  since  two 
grid-stretching  transformations  are  needed  for  the  square  cross  section,  as  opposed 
to  one  for  the  circular  pipe.  It  is  also  worth  noting  that  the  computed  velocity 
contours  in  Fig.  9a  in  the  straight  section  2.5  duct  widths  upstream  of  the  bend 
do  not  display  the  "corner  distortion"  found  in  the  measurements.  This  distortion  is 
the  result  of  weak  "stress  driven"  secondary  flows,  and  its  absense  in  the  computed 
results  is  attributed  to  the  assumption  of  an  isotropic  turbulent  viscosity  in  the 
turbulence  model. 


Flow  Structure  for  Pipe  Bend  Cases 

A  sequence  of  plots  is  given  in  Figs.  10  -  15  for  each  of  the  three  pipe  flow 
cases.  These  figures  are  designed  to  show  the  flow  structure  as  it  develops  at 
successive  axial  locations  for  each  computed  solution,  the  axial  distribution  of 
pressure  coefficient  referred  to  reference  conditions  and  bulk  velocity  at  the  inner 
and  outer  wall  locations  in  the  plane  of  the  centerline  is  first  given.  The  pressure 
coefficient  provides  an  indication  of  the  extent  of  influence  of  the  flow  in  the  bend 
on  the  flow  in  the  upstream  and  downstream  straight  extensions.  The  distribution  of 
grid  points  is  included  in  these  plots.  Also  shown  is  the  difference  in  pressure 
coefficient,  Ac^,  which  would  occur  for  a  potential  flow  with  velocity  inversely 
proportional  to  radius,  in  a  bend  of  the  same  curvature. 

Following  the  pressure  coefficient,  a  sequence  showing  the  primary  and  secondary 
velocity  distributions  at  successive  axial  locations  is  given.  The  same  contour 
values  are  shown  throughout  Figs.  10  -  15;  the  vector  magnitudes  are  renormalized 
for  each  plot  and  thus  indicate  only  flow  direction  and  relative  magnitude  within 
each  plot.  The  strength  of  the  secondary  flows  is  shown  by  giving  the  magnitude  of 
the  peak  velocity,  VMAX,  within  the  vector  plot. 
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Comparing  the  laminar  and  turbulent  cases  in  Figs.  10  -  13,  the  laminar  flow 
is  seen  to  have  a  mush  larger  pressure  drop  (viscous  loss)  and  higher  peak  axial 
velocity  than  the  turbulent  case.  The  secondary  flows  are  of  similar  strength, 
reaching  about  40  per  cent  of  the  mean  axial  velocity.  The  turbulent  pipe  elbow 
case  (Figs.  14,  15)  has  an  inflow  identical  to  that  of  Case  2,  but  there  are 
considerable  differences  in  flow  development,  as  a  result  of  the  strong  curvature. 
The  radial  pressure  gradient  across  the  bend  is  very  large  (Fig.  14),  and  the 
secondary  velocities  are  over  60  per  cent  of  the  mean  axial  velocity.  In  addition, 
the  peak  axial  velocity  is  larger  (  »  1.6)  than  the  2.3  curvature  case,  and  occurs 
near  the  inside  of  the  bend  rather  than  the  outside.  Somewhat  surprisingly,  the 
only  flow  separation  present  was  of  very  limited  extent  and  is  not  visable  in  the 
plots.  This  occurred  near  the  90-degree  location,  on  the  inner  side  of  the  pipe 
bend. 


Influence  of  Boundary  Layer  Thickness  on  Duct  Flow 

The  axial  and  secondary  velocity  distributions  at  60  and  90  degrees  are 
compared  in  Figs.  16  -  17  for  the  two  turbulent  duct  flows  (Cases  5,  6)  with  dif¬ 
ferent  inlet  boundary  layer  thickness,  <5^.  The  computed  flow  structure  is  not  very 
sensitive  to  this  parameter,  and  the  most  notable  difference  is  that  the  distortion 
of  the  axial  velocity  produces  a  thinner  shear  layer  on  the  outer  wall  for  case  6 
with  <5^n  ■=  0.5.  Finally,  the  pressure  coefficient  for  these  two  cases  is  shown  in 
Fig.  18.  The  pressure  field  remains  largely  two  dimensional,  although  there  is  some 
three-dimensional  distortion  associated  with  the  corner  v.rtex  flow  on  the  inside  of 
the  bend. 
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■—  V 


CONCLUDING  REMARKS 


Three-dimensional  laminar  and  turbulent  flow  within  90-degree  bends  of  strong 
curvature  and  both  circular  and  square  cross  section  has  been  studied  by  numerical 
solution  of  the  compressible  Reynolds-averaged  Navier-Stokes  equations.  Straight 
extensions  upstream  and  downstream  of  the  bend  are  included  in  the  computed  flow 
region.  It  is  believed  that  the  physical  processes  involved  require  the  viscous 
sublayer  to  be  resolved,  and  the  computational  approach  provides  for  this 
resolution.  Six  different  flow  cases  were  considered,  and  the  results  were 
evaluated  by  comparison  with  experiment  measurements.  In  general,  very  good  quali¬ 
tative  and  reasonable  quantitative  agreement  between  solution  and  experiment  was 
observed.  The  developing  flow  structure  and  its  dependence  on  geometric  and  flow 
parameters  was  also  studied. 

Collectively,  this  sequence  of  experimental  comparisons  helps  to  establish 
the  accuracy  with  which  these  flows  can  be  predicted  by  the  present  method  using 
moderately  coarse  grids  (  «  10,000  grid  points).  With  some  allowance  made  for  the 
complexity  of  the  flow  and  the  difficulty  of  both  computing  and  measuring  flows  of 
this  type,  the  agreement  between  computed  and  measured  results  is  regarded  as 
generally  very  good.  The  sources  of  disagreement  are  believed  attributable  to 
numerical  truncation  error  and  to  a  lesser  extent  the  turbulence  model. 
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(a)  Plane  of  Centerline 
Curvature 


(b)  Plane  Normal  to 
Centerline 


Fig.  1  -  Schematic  of  Coordinate  System. 
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Fig.  A  -  Convergence  Rate  for  Selected  Cases. 
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Laminar  pipe  Bend  (Case  1  -  Re  «  500,  R/d  «  2.8). 
Axial  Velocity  Compared  with  Measurements  of  [13] 


d  (case  2  -  Re  ! 
elocity  Comparei 


Turbulent  Duct  Bend  (Case  5  -  Re  ■  40,000,  R/d  * 
Axial  Velocity  Compared  with  Measurements  of  [12] 
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Wall  Pressure  Coefficient 
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<E)  90° 


Fig.  11  (Continued) 
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Wall  Pressure  Coefficient 


Grid  Spacing 


Wall  Pressure  Coefficient 


Grid  Spacing 


Axial  Location 

Fig.  14  -  Turbulent  Pipe  Elbow  (Case  3  -  Re  *  43,000, 
R/d  =  1.0).  Wall  Pressure  Coefficient  in 
Centerline  Plane  of  Symmetry. 
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Fig.  15  (Continued) 
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